Reduced-power GPS-based system for tracking multiple objects from a central location

ABSTRACT

Location of an object to be tracked is determined by measuring, at a receiver situated at the object, the propagation time differences between the signals from a plurality of GPS satellites, each of which is received by the receiver situated at the object. The measured propagation time difference values are transmitted to a central station where the location of the object to be tracked is calculated based upon the propagation time differences of the signals transmitted from the satellites and data derived from a receiver apart from the object but also receiving signals from the satellites.

This application is a continuation of application Ser. No. 08/465,229filed May 31, 1995 now abandoned.

FIELD OF THE INVENTION

This invention relates to object tracking systems used to locate andtrack multiple terrestrial (including water-based) objects and, moreparticularly, to power efficient object tracking systems based uponinformation obtained from satellites.

BACKGROUND OF THE INVENTION

The tracking and location of assets such as railcars, shipping or cargocontainers, trucks, truck trailers, automobiles, etc. can be highlyadvantageous in commerce. Precise tracking of such vehicles and objectscan facilitate their being allocated and positioned in an efficientmanner, and can provide for immediate, accurate localization of lost,delayed or damaged assets. The space-based global positioning system(GPS) implemented by the United States Department of Defense constitutesa convenient instrumentality for determining geographical position inreal time.

The GPS is a multiple satellite-based radio positioning system in whicheach satellite transmits data that allows precise measurement of thedistance from selected ones of the GPS satellites to the antenna of auser's receiver so as to enable the user to compute position, velocityand time parameters through known triangulation techniques. The signalsprovided by the GPS can be received both globally and continuously.

The GPS comprises three major segments known as the space, control anduser segments. The space segment consists of 21 operational satellitesand three spare satellites. The satellites are positioned in aconstellation such that typically seven satellites, but a minimum offour, are observable by a user anywhere on or near the earth's surface.Each satellite transmits signals on two frequencies known as L1 (1575.42MHz) and L2 (1227.6 MHz), using spread spectrum techniques that employtwo types of spreading functions.

C/A (or coarse/acquisition code) and P (or precise) pseudo random noise(PRN) codes are transmitted on frequency L1, and P code only istransmitted on frequency L2. The C/A is available to any user, militaryor civilian, but the P code is only available to authorized military andcivilian users. Both P and C/A codes contain data that enable a receiverto determine the range between a satellite and the user.

Superimposed on both the P and C/A codes is a navigation (NAV) message.A NAV message contains the GPS signal transmission time; a handover wordused in connection with the transition from C/A to P code tracking;ephemeris data for the particular satellites being tracked; and almanacdata for all of the satellites in the constellation, includinginformation regarding satellite health, coefficients for the ionosphericdelay model for C/A code users, and coefficients used to calculateuniversal coordinated time (UCT).

The control segment comprises a master control station (MCS) and anumber of monitor stations. The monitor stations passively track all GPSsatellites in view, collecting ranging data and satellite clock datafrom each satellite. This information is passed on to the MCS where thesatellite's future ephemeris and clock drift are predicted. Updatedephemeris and clock data are uploaded to each satellite forretransmission in each satellite's navigation message. The purpose ofthe control segment is to ensure that the information transmitted fromthe satellite is as accurate as possible.

The GPS is intended to be used in a wide variety of applications,including space, air, sea and land vehicle navigation, precisepositioning, time transfer, altitude referencing and surveying. Atypical GPS receiver comprises a number of subsystems, including anantenna assembly, an RF (radio frequency) assembly, and a GPS processorassembly. The antenna assembly receives the L-band GPS signal andamplifies it prior to insertion into the RF assembly. A significantfactor affecting accuracy of the computed position, velocity or timeparameters is the positional geometry of the satellite selected frommeasurement of ranges. Generally, a best position solution is obtainedusing satellites having wide angles of separation. Considerable emphasishas therefore been placed on designing antenna systems to receive, withuniform gain, signals from any point on the hemisphere.

The RF assembly mixes the L-band GPS signal down to a convenient IF(intermediate frequency) signal. Using various known techniques, the PRNcode modulating the L-band signal is tracked through code-correlation atthe receiver. This provides the processing gain needed to achieve asignal-to-noise (SNR) sufficient for demodulating the navigation dataand signal transmission time. The Doppler shift of the received L-bandsignal is also measured through a carrier tracking loop. The codecorrelation and carrier tracking function can be performed using eitheranalog or digital signal processing.

By differencing the signal transmission time with the time of reception,as determined by the clock of the receiver, the pseudo range between thereceiver and the satellite being tracked may be determined. The pseudorange includes both the range to the satellite and the offset of theclock from the GPS master time reference. The pseudo range and Dopplermeasurements (and the navigation data) from four satellites are used tocompute a three dimensional position and velocity fix, which calibratesthe receiver's clock offset and provides an indication of GPS time.

In some known receivers, the receiver processor controller (RPC)functions are performed using a computer separate from that on which thenavigation functions are performed. In other known receivers, both typesof functions are performed by a single computer. The RPC processing andmemory functions performed by a typical GPS receiver include monitoringchannel status and control, signal acquisition and reacquisition, codeand carrier tracking loops, computing pseudo range (PR) and delta range(DR) measurements, determining data edge timing, acquisition and storageof almanac and ephemeris data broadcast by the satellites, processorcontrol and timing, address and command decoding, timed interruptgeneration, interrupt acknowledgment control and GPS timing.

U.S. Pat. No. 5,225,842 describes an apparatus and method for computingthe position and velocity of multiple low cost vehicle-mounted sensors,monitored and tracked by a central control station. The receiverprocessor functions are physically separated from the navigationfunctions and the low rate data interfaces provided between thecomputers that perform these functions, thus achieving cost saving inthe GPS sensor that is employed on board each vehicle.

One type of known GPS receiver is described in U.S. Pat. No. 4,114,155,wherein the position of a receiver responsive to C/A signals derivedfrom multiple, orbiting spacecrafts is determined to an accuracy betterthan 300 meters. Each of the C/A signals has the same carrier frequencyand a different, predetermined Gold code sequence that normally preventsposition determination from being more accurate than to within 300meters. C/A signals transmitted to the receiver are separately detectedby cross-correlating received Gold code sequences with plural locallyderived Gold code sequences. Four of the detected C/A signals arecombined to compute receiver position to an accuracy of 300 meters. Todetermine receiver position to an accuracy better than 300 meters, therelative phase of internally-derived Gold code sequences is varied overthe interval of one chip (i.e., pulse) of each sequence, to derivesecond cross-correlation values indicative of received andlocally-derived Gold code sequences.

The basic approach followed most recently is to receive and process thesignals from several of the GPS satellites in order to determine rangeto each satellite (and relative velocity). With perfect knowledge ofrange to only three of the GPS satellites, exact receiver position canbe determined from the intersection of the three “spheres” induced bythe known satellite positions and the derived receiver ranges. Withreceiver noise and imperfect knowledge of satellite positions, thereceiver-satellite ranges can only be estimated. Typically, errors fromreceiver noise are reduced by (effectively) averaging many rangecalculations.

In the above most recent approach, the range from a particular satelliteis estimated by reading a time stamp from the satellite's data stream(the transmission instant), subtracting this from the reception time,and multiplying the time difference by the speed of light. Any error insatellite and receiver clock synchronization leads to proportional rangeerrors. Because the same clock is used in receiving from all satellites,there is only one unknown receiver clock “bias”. By using a fourth (ormore) satellite, the clock bias and ranges can be jointly estimated.

At the receiver, the reception time is determined by performing across-correlation of the received data with a local replica of the knownsatellite Gold code, and noting the time of a chosen correlation peak,and its position relative to the time stamp. The satellite signalstructures use Code Division Multiple Access (CDMA) so that the abovecross correlation is part of the standard GPS receiver processing.

The above-described system that follows the most recent basic approachassumes that each receiver must determine its own position. In thesystem of the invention, there is a central facility or station thatneeds the receiver positions and can communicate with the receivers.Each tracked object (e.g., a railcar) carries a GPS-based receiver thatprocesses data from several of the visible GPS satellites. However, thefull position determination is not made at the railcar. Instead, onlypartial processing is done at the railcar and intermediate results aretransmitted to the central station. The forms of both the partialprocessing and intermediate results are chosen to minimize thecomplexity and energy requirements at the railcars.

The standard GPS system requires that the transmit-time stamps,satellite ephemeris and other correction data be decoded from eachsatellite's data stream at the tracked object. The receiver is thusrequired to process data from each satellite long enough (between sixand 150 seconds) to synchronize with, and decode, these data. Thisconsumes significant power.

SUMMARY OF THE INVENTION

Briefly, in accordance with a preferred embodiment of the invention, amethod for identifying location of an asset or object to be trackedcomprises measuring data related to propagation time differences betweensignals transmitted from a plurality of satellites and received at theobject to be tracked, transmitting the data to a central station, andcalculating, at the central station, the location of the object to betracked based upon the transmitted data. The data received at the objectto be tracked may include data identifying a respective associatedsatellite by, for example, a satellite identification number, such thatthe step of calculating the location of the object to be tracked isthereupon based further upon the satellite identification numbers.

In another preferred embodiment, a reduced-power GPS-based system fortracking location of an asset or object from a central locationcomprises a central station at the central location, and an object to betracked which includes means for receiving signals from at least fourGPS satellites, first processor means for processing data from thereceiver means propagation time differences for the signals, andtransmission means for transmitting the processed data to the centralstation. Second processor means situated at the central stationdetermine location of the object based upon the data received from thetransmission means.

In still another preferred embodiment, a reduced-power GPS-based systemfor tracking multiple objects from a central location comprises acentral station at the central location, and a plurality of objects tobe tracked. Each of the objects includes receiver means containing anantenna for receiving signals including data related to the propagationtime differences of the signals from at least four GPS satellites, firstprocessor means for calculating a receiver code word phase for each ofthe satellites based upon the signals received by the receiver means,and transmission means for transmitting the code word phase to thecentral station. Second processor means are provided at the centralstation for determining the signal propagation times between theplurality of satellites and each of the tracked objects based on thereceiver code word phase transmitted by the respective object and fordetermining from the signal propagation times the location of eachrespective one of the objects.

In another preferred embodiment, a reduced-power GPS-based system fortracking location of multiple objects from a central location comprisesa central station at the central location, and an object to be trackedwhich includes means for receiving signals including data related topropagation time differences of said signals from at least four GPSsatellites, first processor means for calculating a receiver code-timeoffset for each of the satellite signals and for determining thereceiver code period for each signal, and for determining identificationnumbers of the at least four GPS satellites, and transmission means fortransmitting the receiver code-time offsets, code periods, andidentification numbers to the central station. Second processor means atthe central station determine the signal propagation times between theplurality of satellites and the tracked object and determine location ofthe object based upon the receiver code-time offsets, code period, andsatellite identification numbers transmitted by the transmission means.

In still a further preferred embodiment, a reduced-power GPS-basedsystem for tracking location of an object from a central locationcomprises a central station at the central location, and an object to betracked which includes means for receiving transmission signals from atleast four GPS satellites, first processor means for calculating areceiver bit phase for each of the satellite signals based upon thesignals received by the receiver means, means for keeping track of timeat the object, and transmission means for transmitting the bit phasesignals and time signals to the central station. Second processor meansare provided at the central station for determining the signalpropagation times between the plurality of satellites and the object andfor determining location of the object based upon the bit phase and timesignals transmitted by the transmission means.

Utilizing the present invention, power consumption and calculationcomplexity at the tracked object are reduced relative to that for astandard GPS receiver. Arrival time differences between satellitesignals are measured at the tracked object and this information isrelayed to the central station via the separate communications link.Satellite data streams need not be decoded at the tracked object.

The central station thereby necessarily determines the location of theobject to be tracked. Because the receiver front end and data processoruse significant power only when processing, the receiver power can bedramatically reduced by being “active” only long enough to get accuratetime-difference measurements. This can be less than one second andrequires no GPS data-frame synchronization because of the nature of thesignals. For example, assuming that the tracked object is a railcar, newrailcar locations typically are needed only as frequently as 15 minutes.Thus the receiver energy used is reduced in direct proportion to thereduction of “active” receiver time. Moreover, receiver complexity andcost can be reduced by replacing the advanced microprocessor employed incurrent GPS receivers with a simpler one that is matched to thearrival-time differencing tasks.

In accordance with the invention, one object is is to provide aGPS-based asset tracking system in which processing is performed at alocation remote from the tracked assets and based upon specificrecognizable variables.

Another object is to provide a GPS-based asset tracking system whichrequires minimal energy at the tracked assets.

BRIEF DESCRIPTION OF THE FIGURES

The features of the invention believed to be novel are set forth in theappended claims. The invention, however, together with further objectsand advantages thereof, may best be understood by reference to thefollowing description taken in conjunction with the accompanyingdrawing(s) in which:

FIG. 1 is a block diagram of a remote tracking system in accordance withthe present invention;

FIG. 1A is a block diagram of a railcar tracking unit on an object to betracked in accordance with the invention;

FIG. 1B is a block diagram of the transmitter at the central station ofthe remote tracking system of the invention;

FIG. 2 illustrates the long time-scale data frame relationship between aGPS standard time mark and two satellite data frames in accordance withthe present invention;

FIG. 3 illustrates a typical Gold-code correlation output signal on ashort time scale when the proper code replica is used at the localstation receiver;

FIG. 4 illustrates the time delay relationships of the transmitted andreceived signals;

FIG. 5 illustrates a plurality of intersecting time-difference solutionregions on the surface of the earth;

FIG. 6 illustrates tracked object message timing diagrams in accordancewith the present invention;

FIG. 7 illustrates tracked object location-dependent communicationdelay; and

FIG. 8 is a block diagram of a system for evaluating GPS algorithms inaccordance with the present invention.

DETAILED DESCRIPTION

The present invention is directed to a system and method for reducingthe power and complexity requirements of a local GPS receiver, which canbe carried by a railcar, by effectively requiring only measurement ofarrival-time differences between a plurality of GPS satellite signals.Data related to these time differences are transmitted to a centralstation where the majority of calculations required to determine thereceiver (railcar) location are performed. In the preferred embodiments,a standard CDMA receiver is employed with radio frequency/intermediatefrequency (RF/IF) front end and Gold-code cross correlators.

In FIG. 1, the invention is shown as comprising a plurality of GPSsatellites 12, an object to be tracked, such as a railcar carrying atracking unit 14, and a central station 16. Although the invention isherein described in the context of a railcar, the teachings of theinvention are applicable to a variety of objects which may be tracked bya GPS or satellite-based system. With respect to the GPS signal format,accurate timing of GPS signals is critical and is monitored by centralstation 16.

Each satellite 12 includes its own set of clock correction parameterswithin its data stream. These allow a receiver to ascertain the absolutetransmission timing for each satellite with respect to a common GPSstandard time. A particular satellite's clock may drift relative tothose of other satellites. The GPS system control monitors these offsetsand periodically includes them in the satellite's data stream. The clocktime offsets are not needed at the individual receivers and can bedetermined at central station 16 by utilizing a standard GPS receiverthere or at a remotely controlled site.

As shown in FIG. 1A, railcar tracking unit 14 is actually comprised of arailcar receiver 2 responsive to the signals from GPS satellites 12, aprocessor 3, and a transmitter 4. The received signals are processed inprocessor 3 to ascertain and utilize propagation time differences amongthe signals received from GPS satellites 12. The processed signals arefurnished to transmitter 4 from whence they are transmitted, as reducedorder parameters, to central station 16.

As shown in FIG. 1B, central station 16 includes a receiver 7 responsiveto signals from transmitter 4 of tracking unit 14 (FIG. 1A) and aprocessor 8 responsive to receiver 7 for determining location oftracking unit 14.

In the long time scale data frame shown in FIG. 2, a_(i) is the frametime offset for satellite i. A value a_(j)−a_(i) is the time offsetbetween transmitted signals from two satellites, i and j.

A Gold-code correlator output waveform r_(i)(τ) for satellite i isillustrated in FIG. 3 from the perspective of the railcar receiver. Eachcorrelation-peak location respectively marks the beginning of a new Goldcode cycle in the received waveform. Each Gold code cycle is 1 ms longand comprises 1023 binary Gold code “chips”. Furthermore, there are 20code cycles for each encoded navigation data bit. FIG. 3 alsoillustrates, by a dotted line, a typical bit-boundary position.

With respect to a particular time t_(R) at a railcar receiver, thereceiver code-time offset for satellite i is γ_(i). The receivercode-time offset is the time elapsed to a time t_(R) from the beginningof the code word (cycle) in which t_(R) falls. Similarly, the receiverbit-time offset β_(i) is the time elapsed to time t_(R) from thebeginning of the bit in which t_(R) falls. The satellite-railcar radialvelocity component varies for different satellites and this results in arelative waveform expansion or compression (Doppler) at the railcar.Thus the observed code and bit periods are satellite dependent. The codeand bit periods observed at the railcar for satellite i are designatedas T_(i) ^(C) and T_(i) ^(B), respectively.

Frequently, the railcar receiver will use satellite signals that are notvisible from (i.e., cannot be received by) the central station. Thispresents no problem because the satellite clocks drift slowly (less thanfive nanoseconds error per hour). If, at the railcar receiver, timedifferences are calculated using a satellite that is not visible at thecentral station, then, at the central station, the last calculated clockoffset for that satellite can be used (or extrapolated, based on pastdrift rate) until that satellite is again visible. As an alternative,central station 16 (FIG. 1) could communicate with standard GPSreceivers strategically situated to guarantee satellite visibility.

A key feature of the present invention is the provision of a method fordetermining location of the object (here, a railcar) to be tracked. In afirst method (“method 1”), the object's location is accuratelydetermined from propagation time differences between at least fivesatellite' signals received at the tracked object. This method requiresno measurement of time at the tracked object. The propagation timedifference between signals from satellites i and j is defined asΔ_(ij)=τ_(j)−τ_(i), where τ_(i) is the signal propagation time fromsatellite i to the railcar. These propagation times are not directlymeasured at the receiver but are calculated from received code word orbit phases, as described below (see equations 8 and 9). Thesatellite-railcar ranges and propagation time differences are related bythe equationR_(i)(t, t−τ_(i))−R₁(t, t−τ₁)=CΔ_(1i)(t),   (1)or:R₂(t, t−τ₂)−R₁(t, t−τ₁)=CΔ₁₂(t)R₃(t, t−τ₃)−R₁(t, t−τ₁)=CΔ₁₃(t)R₄(t, t−τ₄)−R₁(t, t−τ₁)=CΔ₁₄(t)R₅(t, t−τ₅)−R₁(t, t−τ₁)=CΔ₁₅(t)where R_(i)( ) is a function of the parameters in its argument and isdiscussed below. In equation (1), t is the common time at which thesignals are measured at the receiver, C is the speed of light, andR ₁(t,t−τ ₁)=[(x(t)−x _(i)(t−96 _(i)))²+(y(t)−y _(i)(t−τ _(i)))²+(z(t)−z_(i)(t−τ _(i)))²]^(1/2)   (2)is the range from the i'th satellite to the railcar receiver. Also,x_(i), y_(i) and z_(i) are time-dependent coordinates for satellite iand are specified by the satellite's ephemeris equations. For satellitei, the signal received at time t=t_(R) was transmitted at timet_(R)−τ_(i), which is now defined as t_(i) ^(T). Time t_(R) is identicalfor all of the satellites. Furthermore, the propagation times arerelated byτ_(i)=τ₁+Δ₁₁.   (3)

For any particular measurement time t=t_(R), equation (1) can berewritten asR ₂(x,y,z,t ₁ ^(T)−Δ₁₂)−R ₁(x,y,z,t ₁ ^(T))=CΔ ₁₂   (4)R ₃(x,y,z,t ₁ ^(T)−Δ₁₃)−R ₁(x,y,z,t ₁ ^(T))=CΔ ₁₃R ₄(x,y,z,t ₁ ^(T)−Δ₁₄)−R ₁(x,y,z,t ₁ ^(T))=CΔ ₁₄R ₅(x,y,z,t ₁ ^(T)−Δ₁₅)−R ₁(x,y,z,t ₁ ^(T))=CΔ ₁₅The corresponding satellite transmission times are needed to properlydetermine the satellite locations needed in equation (2). From equations(2), (3) and (4), the object coordinates x, y, z and the transmissiontime t₁ ^(T) form the only unknowns. Yet the time delay differences inequation (1) do not have to be calculated with reference to satellite(1); i.e., any satellite pairing will be satisfactory as long as thefour equations in (1) utilize unique pairings.

Using these simultaneous equations at central station 16 (FIG. 1), fourtime delay differences determine a point at the intersection of fourhyperbolic sheets in three dimensional space as well as the unknowntransmit time, t₁ ^(T). The other satellite transmission times t₂ ^(T),t₃ ^(T), t₄ ^(T), and t₅ ^(T) are then found from the equation:t ₁ ^(T) =t ₁ ^(T) −Δ1i.   (5)

The resultant values are utilized to determine the satellite locationsneeded in equation (2). The nonlinear equations (4) are readily solvedwith standard techniques using iteration. It is not necessary that timet_(R) be measured at the railcar receiver. This value can be found, ifdesired, at the central station from the equation:t _(R) =t ₁ ^(T) +R ₁ /C.   (6)

In the present invention, the propagation time differences Δ_(ij) arenot directly measured at the railcar receiver. Instead, only the code orbit phases associated with reception time t_(R) are measured, and thesedata, or their differences including the satellite identifications, aretransmitted to the central station. This permits the railcar receiver tofocus only on cross correlation and to operate long enough to get asufficient SNR (Signal-to-Noise Ratio) through averaging. In general,the averaging times are so brief that any Doppler shifts can beconsidered as constants. Thus the bit and code periods are constantduring these periods. Railcar locations can be determined at the centralstation by calculating the true differential delays from the measuredphases, as described below, and using them in equations (1) through (5).

In determining the propagation time difference between two satellitesignals, it is assumed that the Doppler shift is constant locally, i.e.,that the satellite-railcar propagation delay changes linearly with time.Using this assumption, a particular relative position in a received codeword or data-bit cycle corresponds to the same relative position in theassociated transmitted code word or data-bit cycle. From the previousdefinitions, and as shown in FIG. 4, for each satellite i, j, etc.,t _(R)=t₁ ^(s) t+T ^(C)(m _(i) +μ _(i))+τ_(i).   (7)Pursuant to equation (7), t_(i) ^(s) is the time when the start of acurrently received frame was transmitted, m_(i) is the integer number ofcode periods between time t_(R) and the beginning of the received dataframe, and μ_(i) is the receiver code word phase at time t_(R) and isdefined asμ_(i)=γ_(i) /T _(i) ^(C)   (8)where γ_(i) is the receiver code-time offset at the receiver, and T_(i)^(C) is the code period (at time t_(R)) in the received data frame.T^(C) is the common code period at all transmitters. Because the Dopplershift is constant, it is true thatε_(i) /T ^(C)=γ_(i) /T _(i) ^(C)   (9)so the relative code word phase in the transmitted waveform at timet_(R)−τ_(i) is equal to μ_(i) that was measured at the receiver at timet_(R). From equation (7), the propagation time difference betweensignals from satellites i and j isΔ_(ij)=τ_(j)−τ_(i) =a _(j) −a _(i) +T ^(C)(m ^(i) −m _(j))+T^(C)(μ_(i)−μ_(j))   (10)which utilizes the equation t_(i) ^(S)=t_(GPS)+a_(i) for each i. At thecentral station, a_(i), a_(j) and T_(C) are known, t_(GPS) is the chosenreference time, and μ_(i) and μ_(j) are received from the tracked objectas measurements. The integers m_(i) and m_(j) are unknown and yieldambiguity as discussed below. At the receiver, time t_(R) can be chosento align with the first satellite's received code word boundary. Thenμ₁=0 by convention and the other four phases can be sent to the centralstation. In this way, the central station can know the individual phasevalues as well as their differences. Knowledge of each phase value atthe central station constrains each unknown associated transmission timeto be on a lattice of time points that have that phase value.

While the above discussion illustrates how the several satellitepropagation time differences can be deduced from receiver code wordphase measurements, as an alternative, receiver bit phases can bemeasured. In this instance, μ_(i) of equation (8) becomes the bit phasewhen γ_(i) is replaced by β_(i), and T_(i) ^(C) is replaced by T_(i)^(B) such that μ_(i)=β_(i)/T_(i) ^(B). A relationship similar toequation (9) also holds for bit phases. Finally, in equation (10), m₁and m_(i) become the unknown integer number of bit periods. In eitherinstance, the unknowns m_(i) and m_(j) cause a periodic ambiguity inpropagation time difference Δ that must be resolved. This ambiguity hasa period of approximately one millisecond for code periods and 20milliseconds for bit periods.

Because each possible propagation time difference value induces athree-dimensional hyperbolic sheet for the railcar location solution,the ambiguity in Δ induces multiple sheets for each value of Δ. On anassumed flat earth, the sheets from just one phase difference form a setof hyperbolas with the positions of satellites i and j as the foci.Phase differences from other satellite pairs induce difference hyperbolasets. The only feasible railcar location solutions are those where thereis a congruence of one hyperbolic sheet from each participatingsatellite pair. With ambiguous time differences from four satellites,only a small set of possible joint ambiguity images yield intersectinghyperbolas near the earth's surface, and the set is further reduced asmore satellites are used.

Time-delay difference and satellite location accuracies are such thatthe location solution regions from each satellite pair are very narrowrelative to the ambiguity spacing, as shown in FIG. 5. Thus the severalambiguous time delay differences yield a small set of possible railcarlocations.

The use of bit-time differences provides advantages over code-timedifferences because the former yield a smaller number of possiblelocation solutions. The 20 millisecond period for bit-time differencesalways yields a location ambiguity spatial period of at least 1500 milesat the earth surface, while for code-time differences the ambiguityspatial period is only 75 miles. This shortest ambiguity periodassumption is derived by assuming the two satellites involved are on thehorizon in opposite directions. For railcar tracking applications, carlocations will be known a-priori to within 1500 miles, so a 1500 mileambiguity is not a problem. Receiver bit phases, however, are slightlycomplicated to derive because the satellite's binary data stream is NRZcoded.

Using an NRZ code sequence, a sequence of identical bits has nottransitions and bit locations are not visible. However, with GPS data,such sequences are very short, so that bit edges are readily observed.Once a single bit transition has been observed, the receiver bit phasefor the chosen time t_(R) is readily deduced because the code period isobserved and there are always 20 code cycles per bit period.Furthermore, each bit boundary coincides with a code word boundary.

In order to solve equation (4) for the railcar receiver location, thepropagation time delay differences are first found, as from equation(10). In equation (10), (m_(i)−m_(j)) is the unknown integer part of thecode or bit period offset between the signals received from satellites iand j. For a given measurement of μ_(i) and μ_(j), equations (10) and(4) yield a different location solution for each value of (m_(i)−m_(j)).Conceptually, each integer value of (m_(i)−m_(j)) must be tried, and theresulting position solution must be checked against known bounds forvalidity. In earthbound (e.g., railcar) applications, altitude is asimple bound against which each position solution can be checked. Also,railcar velocities are constrained, so that new locations cannot beextremely different from previous locations.

Most of the ambiguous solutions are at invalid altitudes, far from theearth's surface. To avoid wasteful calculation, it is desirable todirectly limit the values of (m_(i)−m_(j)) that are used in equation(10). The GPS satellite-earth geometry constrains each time t_(i) to bebetween 58.5 and 79.9 milli-seconds (ms). This is true because GPSsatellites orbit at approximately 25×10³ km above the earth's center sothat the delay from directly overhead is 58.5 ms while the delay from asatellite on the horizon (nearly 4000 miles farther away) is 79.9 ms.The delay difference between two satellites is therefore constrained tothe interval [−21.4, 21.4] milliseconds. At the central station,equation (10) is readily used to deduce values of (m_(i)−m_(j)) thatsatisfy this interval constraint. If code phases are measured, there areapproximately 43 feasible values of (m_(i)−m_(j)) for each measured(μ_(i)−μ_(j)). Thus there are 43 values for each Δ_(ij) in equation (4).If bit phases are measured, there are only 2 or 3 feasible values foreach propagation time difference Δ_(ij) and the resulting list ofpossible railcar locations will be much shorter.

To summarize the foregoing object-tracking method (“method 1”), fivesatellite signals must be received. Five receiver code word or bitphases are measured, as are the associated satellite identificationnumbers, and these are sent to the central station. From thesemeasurements a list of feasible railcar locations is determined at thecentral station. The basic steps are:

1. Phases μ₁ through μ₅ are measured along with their correspondingsatellite numbers and these data are sent to the central station. Themeasured phases can be code word phases in the simplest receiver, ordata-bit phases in a slightly more complex receiver. Bit phase could bespecified as code word phase plus an integer number of code words offsetfrom the bit transition.

2. A standard GPS receiver at the central station determines satellitetransmission offsets a₁ through a₅ that are valid around the time thereceiver measurements are received at the center. Validity of eachsatellite ephemeris equation is determined for the same time period.

3. At the central station, reasonable values for the integer offsets(m₁−m_(j)) for j from 2 through 5 are selected, and equation (10) isused to calculate Δ_(1j) for j from 2 through 5.

4. At the central station, an initial value for t₁ ^(T) is chosen andequation (5) is used to find corresponding values for t_(i) ^(T) for ifrom 2 through 5. The value of t₁ ^(T) is constrained such that thetransmitted bit or codeword phase at t₁ ^(T) has the measured phase, μ₁.The integer m₁ is indirectly specified when choosing a t₁ ^(T) thatsatisfies the above constraint. At the central station, the initialvalue for transmit time t₁ ^(T) may be set to an approximation of thereceiver message time t_(R), if known. In a railcar-trackingapplication, time t_(R) could be known to within several minutes at thecentral station (without any communication of time values).

5. Standard iterative methods are used with equations (1) and (2) inorder to solve for the railcar receiver (x, y, z) position andtransmission time t₁ ^(T). In equation (2), the ephemeris equation foreach satellite is included in the iteration. The last reported locationcan be used as the initial values for x, y, and z in the iteration.

6. Steps 3 through 5 are then repeated for each feasible combination ofinteger offsets (m₁−m_(j)) from 2 through 5. This yields a list ofpotential railcar receiver location solutions.

An advantage of the tracking method described above is that requirementsfor receiver clock accuracy are minimal. Railcar time is not part of themeasurement set.

The present invention may alternatively employ a second method (“method2”) for determining location of an object being tracked. The secondmethod is similar to the first method, except that only four satellitesare used and the railcar receiver message time t_(R) must be made knownto the central station, (e.g., measured and transmitted to the centralstation). By using four satellites, only three independent propagationtime differences can be obtained from the receiver code word orbit-phase measurements. To determine the railcar location, the time atwhich these propagation time differences are valid must be known. Morespecifically, as shown in FIG. 4, the transmission times t_(i) ^(T), fori from 1 to 4, must be determined so that the satellite position(associated with the common reception time t_(R)) can be determined. Ina manner similar to the first tracking method, the railcar receiverlocation is determined from the following three equations:R ₂(x,y,z,t ₁ ^(T)−Δ₁₂)−R ₁(x,y,z,t ₁ ^(T))=CΔ ₁₂   (11)R ₃(x,y,z,t ₁ ^(T)−Δ₁₃)−R ₁(x,y,z,t ₁ ^(T))=CΔ ₁₃R ₄(x,y,z,t ₁ ^(T)−Δ₁₄)−R ₁(x,y,z,t ₁ ^(T))=CΔ ₁₄.In equation 11, the Δ_(1j) terms are derived from the measured code wordor bit phases with ambiguities as in the first tracking method. Also asin the first method, (t₁ ^(T)−Δ_(1i)) is the transmission time from thei'th satellite corresponding to reception time t_(R). This transmissiontime determines the satellite location from the ephemeris equations, andthe satellite locations are needed to determine the ranges in equations(11). The Δ_(1i) values are known at the central station from thereceived measurements and equation (10). Furthermore, time t₁ ^(T) isrelated to time t_(R) byt ₁ ^(T) =t _(R)−τ₁(x,y,z,x ₁(t ₁ ^(T)),y ₁(t ₁ ^(T)),z ₁(t ₁ ^(T))).  (12)Here τ₁ is the GPS signal propagation delay from satellite 1 anddepends, as shown, on location of the object being tracked, andsatellite location. The satellite location (x₁, y₁, z₁) depends, inturn, on time t₁ ^(T). Therefore, if time t_(R) and the satelliteephemeris equations are known, then time t₁ ^(T) depends only on thelocation of the object being tracked. For a given location of the objectbeing tracked, equation (12) can be solved iteratively for t₁ ^(T).After convergence, the value of t₁ ^(T) can be modified to the nearestpoint of the time grid induced by equations (8) and (9) and the known(measured) value of μ₁ so as to obey the previously mentioned latticeconstraint. To speed convergence, this constraint can be applied aftereach iteration. Equations (11) and (12) together form a system ofnonlinear equations with object location (x, y, z) as the only unknowns,and these equations can be solved using standard iterative techniques.

In many railcar tracking applications (e.g., when the railcar istraveling), extreme accuracy is not required. In such cases, t_(i) ^(T)can be approximated by t_(R)−(79.9−58.5)/2. As stated earlier in thedescription of the first tracking method, the GPS-signal transmissiontimes (t_(i) ^(T)) must lie between t_(R)−79.9 ms and t_(R)−58.5 ms.This yields at t₁ ^(T) uncertainty of 21.4 ms if t_(R) is known. Thesatellites travel at approximately 3.49 km per second so that a t_(i)^(T) uncertainty of 21.4 ms yields a satellite location (x, y or z)uncertainty of, at most, 74.8 meters, which translates directly to asimilar uncertainty in railcar receiver location. This may providesufficient accuracy for most railcar-tracking applications. If greateraccuracy is desired, equation (12) can be included in the iterativesolution of equation (11) to determine location of the object beingtracked. The initial value for t₁ ^(T) can be set equal to t_(R) andequation (11) iterated to solution for (x, y, z). This result and theinitial value for t₁ ^(T) can then be used to initialize an iterationfor solving equation (12). This yields an improved value for t₁ ^(T)which can be used with the previous (x, y, z) solution to initializeanother iterative solution of equation (11). This cascade of iterationswill rapidly converge to the correct railcar location, subject to the(m_(i)−m_(j)) induced ambiguities.

Small errors in t_(R) can be tolerated because t₁ ^(T) must obey thepreviously-mentioned lattice constraint. Thus if the received codeworkor bit phases μ_(i) are measured at time t_(R), but the measured time isactually reported as t_(R) plus some error ε, the railcar location willbe accurately calculated for time t_(R) if the value of ε is less thanone-half of the code word (data-bit) period. If desired, the correctvalue for t_(R) can be determined following completion of the iterationusing the following equation:t _(R) =t ₁ ^(T) +T ₁ /C.   (13)

In performing this method of tracking an object, it is assumed that thereception time t_(R) for the railcar receiver can be accurately measuredand used to find the relevant satellite transmission times. This iseasily done if the object being tracked includes an accurate clock andits time reading at t=t_(R) is sent to the central station.Alternatively, time t_(R) can be determined at the central station bydetermining the time of data arrival at the central station and usingthis as a reference. As shown in FIG. 6, data transmitted from a GPSsatellite at time τ_(O) is detected at the railcar receiver at timeτ₁=t_(R), is transmitted to the central station at time τ₂, and isreceived by the central station at time τ₃. The central stationmaintains an accurate clock so that the value t_(R) can be determined ifboth (τ₃−τ₂) and (τ₂−τ₁) are known. The delay (τ₂−τ₁) can be easilyrecorded by an inexpensive timer located at the railcar receiver and canbe transmitted as data to the central station. The communication delay(τ₃−τ₂) is dependent upon the data communication network and may beknown at the central station where reception time t_(R) may then becalculated using the formula:t _(R) =τ ₃−(τ₃−τ₂)−(τ₂−τ₁).   (14)

As indicated in FIG. 7, there are some situations in which the datacommunication may be made via a satellite link. Assuming the use of ageostationary satellite 19 for the data link, transmissions out to andback from the satellite will experience a delay which depends on thelocations of central station 16, railcar receiver 14 and satellite 19.This delay will be bounded between φ=2R₁/C and φ₂=2(R₁+R₂)/C as can beseen from FIG. 7, where C is the speed of light, R₁ is the distance fromsatellite 19 to the surface of the earth and R₂ is the earth's radius.With R₁ approximately 38.623 km and R₂ approximately 6,437 km, the delayuncertainty is approximately 46 ms. During this interval, the GPSsatellites move 174 meters at their speed of 3.78 km per second. Thusthe railcar locations will have an additional error if this delayuncertainty is not resolved.

If a more accurate railcar location is required, the uncertainty ingeostationary satellite delay (due to uncertain railcar location) can beeliminated with an iteration scheme such as that discussed earlier (forthe GPS signal delays). Initially, the railcar location is calculated atthe central station by assuming a particular feasible value for thecommunication delay (τ₃−τ₂). With the new railcar location and the knownsatellite 19 and central station 16 positions, the communication delayis recalculated. This corrected delay is then used in equation (14) tocorrect the value t_(R). The new value of t_(R) is used, as describedearlier, to find a railcar location. This iteration can continue untillittle change in railcar location is observable.

For accurate results from this approach, the data storage delay (τ₂−τ₁)must be short enough such that the railcar receiver clock yields anaccurate storage delay measurement. In a preferred embodiment, therailcar will collect time-difference data approximately every 15 minutesand report these data to the central station hourly. Thus the storagedelay (τ₂−τ₁) for the first measurement set could be up to 60 minutes,or 3600 seconds. An inexpensive receiver clock may keep time to one partin one million, so that a 3.6 ms error may accrue in 3600 seconds. Thiserror is directly reflected in the measured value of t_(R). If theresidual error in t_(R), due to receiver clock drift, is kept to lessthan one-half of the code word or data-bit period, the iteration methoddescribed above will reduce the residual railcar location error to zero.

To summarize this second object-tracking method (“method 2”), foursatellite signals must be received. Four receiver code word or bitphases are measured, as are the associated satellite identificationnumbers and the single measurement time, and these are sent to thecentral station. The basic steps are:

1. Code word or bit phases μ₁ through μ₄ and their associated satellitenumbers are measured at reception time t_(R). These data, including themeasurement time t_(R), are sent to the central station. The bit phasecan be specified as code word phase plus an integer number of code wordsoffset from a bit transition. If an accurate clock is not available atthe railcar receiver, the data-storage time (τ₂−τ₁) at the railcarreceiver is sent instead of the actual measurement time.

2. If storage time is measured at the object being tracked and sent tothe central station, then the time of data reception (τ₃) at the centralstation is determined and used with storage delay (τ₂−τ₁), and with theknown bounds on transmission delay between the object being tracked andcentral station (τ₃−τ₂), to estimate the measured common reception timet_(R) from equation (14).

3. A standard GPS receiver is used at the central station to determinesatellite transmission offsets a₁ through a₄ and the satellite ephemerisequations that are valid near the measured reception time t_(R). Becausethese parameters change relatively slowly, and because the delay betweena GPS satellite and the object being tracked is short, these parameterswill be valid at the several satellite transmission times.

4. At the central station, reasonable values for the integer offsets,(m₁−m_(j)) for j from 2 through 4 are selected, in view of limitationson these offsets that are induced by satellite-earth geometry, priorobject location, object velocity limits, etc. Equation (10) is used tocalculate propagation time difference Δ_(1j) for j from 2 through 4.

5. At the central station, an initial value for time t₁ ^(T) is chosenbased on μ₁, and on limitations on the communication delay (τ₃−τ₂) andbounds on the GPS signal travel delays (τ_(i)). An effective simpleinitial value is t₁ ^(T)=t_(R).

6. Equation (5) is then used to calculate corresponding values for t_(i)^(T) for i from 2 through 4.

7. Standard iterative methods are then used to solve equation (11) forrailcar receiver location (x,y,z). Because the values t_(i) ^(T) arefixed, the ephemeris equations are evaluated only once during thisiteration.

8. For higher accuracy, the new (x, y, z) result is used with equation(12) to iteratively solve for an improved value of t₁ ^(T). This can bedone either after step 7 has converged, or after each iteration of step7.

9. Steps 6 through 8 are repeated, terminating after step 7 if the newposition (x, y, z) is substantially unchanged.

10. Steps 4 through 9 are repeated for each feasible combination ofinteger offsets, (m₁−m_(j)) for j from 2 through 4. This yields a listof potential railcar location solutions.

An advantage of this second method of object tracking is that only foursatellite signals must be received. Furthermore, the iteration forfinding railcar location is less time-consuming because the search forthe correct GPS transmission times is substantially eliminated.

For some communication methods, the communication delay time cannot bedetermined with sufficient accuracy. In such cases, other methods fordetermining receiver and transmission times must be used. One approachis to utilize a unique mark (unique over a sufficient period) in any oneGPS satellite's data stream as a time reference at the tracked object,and to measure reception or receiver time t_(R) relative to thisreceived mark. Because the mark is unique, it can be found in the GPSsignal received at the central station and its transmission time can bedetermined accordingly. Hence the transmission time associated withreception time t_(R) can be found by adding the measured offset to themark time. All other satellite transmission times t_(i) ^(T) can then befound from the known offsets a_(i) at the central station. Utilizingequations (11) and (13), the railcar receiver locations and associatedreceiver times are calculated.

The GPS telemetry-word preamble (TWP) is a specific eight-bit sequence(10001011) that is transmitted at the beginning of every six-secondsubframe from each satellite and is a standard synchronization mark.This sequence cannot be falsely replicated by any prefix or postfix ofup to six bits in length. If a short periodic receiver window issynchronized with a signal event such as the TWP from one of thesatellites whose data is to be processed, and this window is used tooccasionally “awaken” or activate the railcar receiver forTWP-synchronized GPS signal processing, little or no extra power will berequired at the railcar receiver to support learning the transmissiontimes at the central station. This TWP windowing scheme is useful if theuncertainty in satellite transmission times is less than six seconds.Initial window synchronization may require up to a full six seconds ofGPS data processing to ensure TWP acquisition but, once acquired, thewindow can easily be tracked by noting the position of each TWP in itssurrounding window and altering the window timing to keep the TWPcentered. With inexpensive clocks at the railcar receiver, the windowcan drift on the order of only 3.6 ms per hour. If the signal TWP regionis also used for the time-difference processing, then the window candrift only on the order of 3.6/4 ms over the 15 minute inter-measurementinterval. Since the TWP sequence is unique within ±6 bit periods (120ms), use of a 100 ms window ensures that the TWP will not be missed orfalsely recognized.

If the transmission delay to the central station has uncertainty greaterthan six seconds, the TWP windowing scheme can be augmented with GPStime-stamp decoding at the receiver. The GPS time stamp is encoded inthe data stream at a fixed short delay from the TWP word. The TWP windowcan be used to enable time-stamp decoding, and the time stamp can besent to the central station as an estimate of the receiver time t_(R).As an alternative, time t_(R) can be chosen at the railcar receiver tocoincide with the TWP word boundary on satellite channel number 1, andthe following time stamp can be decoded and sent to the central station.The correct value for time t₁ ^(T) is then easily found at the centralstation as the value of the time stamp prior to that from the receiver.

The above railcar receiver windowing scheme need be implemented only onone of the satellite signals in order to determine receiver time t_(R).However, if a separate TWP receiver window is formed, respectively, foreach satellite signal, then all time-delay ambiguities associated withcode or bit periods are readily resolved by noting the TWP-relative timeat the end of the correlation process for each signal. This simplifiesprocessing at the central station at the expense of extra windowprocessing and receiver “on” time at the railcar receiver.

A further method for determining time is to broadcast time signals via aseparate channel accessible to the railcar receivers. For example, timesignals may be transmitted over a separate geostationary satellite linkon a one-second grid or smaller, thereby easing railcar receiver clockaccuracy requirements. This method results in a railcar-locationdependent delay from the central station to the railcar receiver. Thedelay can be accounted for by using iterations in a manner similar tothat described above.

FIG. 8 illustrates a system for evaluating GPS-based localizationalgorithms for low power tracking applications. The system includes astandard GPS hardware unit 22 for receiving GPS signals and developingstandard GPS navigation solutions for reference, and a data interfaceunit 24 and computer workstation 26 for generating GPS navigationsolutions using the reduced-power GPS methods described herein. StandardGPS hardware unit 22 uses a commercially available multi-channel GPSreceiver and GPS signal processor to develop standard navigationsolutions. Data interface unit 24 collects intermediate GPS signalinformation from hardware unit 22 (e.g., demodulated but unprocessed GPSdata streams), performs some data processing (e.g., calculate relativebit phases of the different GPS signals), and passes the results tocomputer workstation 26. The data interface unit also collects decodedGPS almanac, ephemeris, and clock correction data from the hardware unitand passes these to the workstation. At the workstation, the variousreduced-power GPS tracking algorithms are developed using thecommercially available Matlab programming language. Matlab is availablefrom The Math Works, Inc., Natick, Mass. The tracked unit (railcar) andcentral-station functions are performed separately in the computerworkstation.

While only certain preferred features of the invention have beenillustrated and described, many modifications and changes will occur tothose skilled in the art. It is, therefore, to be understood that theappended claims are intended to cover all such modifications and changesas fall within the true spirit of the invention.

Attached hereto, as Appendix A, is a source code listing which may beutilized by the central station to calculate the location of the objectbeing tracked, pursuant to method 1 herein.

APPENDIX A Source Code Listing % ROGPS - Solves for user position andtime based on ephemeris equations % and on time-difference measurementsfor 5 GPS satellites. % U (3×1) the users position % time (1×1) seconds% Ub (3×1) an initial guess at the users position % dU (3×1) is theadjustment made on the initial guess % sat (3×5) is the location of eachsatellite % H (4×4) matrix differentiating r_diff wrt (U,t) % sat0 isthe (3×5) matrix of 5 sat positions at beginning of % interpolationinterval. % sat2 is the (3×5) matrix of 5 sat positions at end ofinterpolation interval. % V is a (3×5) matrix of velocities for all 5satellites. dU = 0; % step 1: Guess an initial Ubar, and time. U_true =[1300748.12;−4500694.5;4313770.5]; % Actual user location Ub =[1200748.12;−4400694.5;4319770.5]; % User location guess time_true =483738.75; time = 483638.75; interval = 300; % Linear interpolation(half) interval for sats. % select five satellites % Files that containephemeris constants. These happen to be % visible at the data collectiontime. sv1 = sv7_eph; sv2 = sv2_eph; sv3 = sv9_eph; sv4 = sv4_eph; sv5 =sv12_eph; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Synthesizemeasurement vector % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % For now,ignore that satellites move during signal propagation % time. Thus, setephemeris time to measurement time, time_true, sv1(1) = time_true;sv2(1) = time_true; sv3(1) = time_true; sv4(1) = time_true; sv5(1) =time_true; % Get satellite positions sat(:,1) =satellite_ephemeris(sv1); sat(:,2) = satellite_ephemeris(sv2); sat(:,3)= satellite_ephemeris(sv3); sat(:,4) = satellite_ephemeris(sv4);sat(:,5) = satellite_ephemeris(sv5); % Calculate perfect measured rangedifferences. for i = 1:4  r_diff_t(i) = norm(sat(:,1) − U_true) −norm(sat(:,i+1) − U_true); end % Add noise to the measurementslight_speed = 3e8; % Meters/sec dt_sigma = 1e-7; % One tenth of a chip.% r_diff_m = r_diff_i + (light_speed * dt_sigma * randn(1,4)); r_diff_m= r_diff_t + [0 0 0 10]; r_diff_m = r_diff_m′; dr_var = dt_sigma * 3e8 *dt_sigma * 3e8; r_diff_covar = [ dr_var 0 0 0; . . . 0 dr_var 0 0; . . .0 0 dr_var 0; . . . 0 0 0 dr_var] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Find satellite average velocities %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Set time to beginning ofinterpolation interval % First element of svi is always time. sv1(1) =time − interval; sv2(1) = time − interval; sv3(1) = time − interval;sv4(1) = time − interval; sv5(1) = time − interval; % Get satellitepositions satA(:,1) = satellite_ephemeris(sv1); satA(:,2) =satellite_ephemeris(sv2); satA(:,3) = satellite_ephemeris(sv3);satA(:,4) = satellite_ephemeris(sv4); satA(:,5) =satellite_ephemeris(sv5); % Set time to end of interpolation interval.sv1(1) = time + interval; sv2(1) = time + interval; sv3(1) = time +interval; sv4(1) = time + interval; sv5(1) = time + interval; % Getsatellite positions. satB(:,1) = satellite_ephemeris(sv1); satB(:,2) =satellite_ephemeris(sv2); satB(:,3) = satellite_ephemeris(sv3);satB(:,4) = satellite_ephemeris(sv4); satB(:,5) =satellite_ephemeris(sv5); % calculate x,y,z velocities (in m/s) V =(satB−satA)/interval; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Forminitial H matrix of partials % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Partials taken at guessed time and user position. % Need partials of satcoordinates w.r.t. time. Over the “interval” % specified above we assumesat velocity is constant. Set partials % to sat velocities. Computevector and range from object to each % satellite. % Get satellitepositions at guessed time. sv1(1) = time; sv2(1) = time; sv3(1) = time;sv4(1) = time; sv5(1) = time; sat(:,1) = satellite_ephemeris(sv1);sat(:,2) = satellite_ephemeris(sv2); sat(:,3) =satellite_ephemeris(sv3); sat(:,4) = satellite_ephemeris(sv4); sat(:,5)= satellite_ephemeris(sv5); % Compute estimated sat user vectors andranges. for i = 1:5, sat_U(:,1) = sat(:,i) − Ub; range(i) =norm(sat_U(:,i)); DV(:,i) = sat_U(:,i)/range(i); end % Form H matrix ofpartials at estimated user location and time; % i'th row of H is H =(pD1i/px, pD1i/py, pD1i/pz, pD1i/pt) % First three elements is just thedifference between the normalized % direction vectors to the first andto the i'th satellites. % Last element is difference between 1st andi'th inner products % of direction vectors with satellite velocityvectors. for i = 1:4, H(i,1:3) = DV(:,i+1)′ − DV(:,1)′; H(i,4) =DV(:,1)′ * V(:,1) − DV(:,i+1)′ * V(:,i+1); end % now compute the firstestimate of G % G = inv(H); loop = 0; epsilon = 1; %meter while (loop <20) loop = loop + 1; % step 3 % Substitute t into ephemeris data and getpositions at guessed time. % Here we use actual ephemeris equations tofind sat positions, % but note that the time derivatives of sat positionare taken from % the linear interpolation. sv1(1) = time; sv2(1) = time;sv3(1) = time; sv4(1) = time; sv5(1) = time; sat(:,1) =satellite_ephemeris(sv1); sat(:,2) = satellite_ephemeris(sv2); sat(:,3)= satellite_ephemeris(sv3); sat(:,4) = satellite_ephemeris(sv4);sat(:,5) = satellite_ephemeris(sv5); sat1(:,loop) = sat(:,1);sat2(:,loop) = sat(:,2); sat3(:,loop) = sat(:,3); sat4(:,loop) =sat(:,4); sat5(:,loop) = sat(:,5); % Compute vector and range fromobject to each satellite. % Needed for new H matrix. for i = 1:5,sat_U(:,i) = (:,i) - Ub; range(i) = norm(sat_U(:,i)); DV(:,i) =sat_U(:,i)/range(i); end % compute r_diff using Ubar and t r_diff(1) =range(1) − range(2); r_diff(2) = range(1) − range(3); r_diff(3) =range(1) − range(4); r_diff(4) = range(1) − range(5); % step 4 % computerange-difference error. r_diff_error(:,loop) = r_diff_m − r_diff′; %step 5 % Form H matrix of partials; % i'th row of H is H = (pD1i/px,pD1i/py, pD1i/pz, pD1i/pt) % First three elements is just the differencebetween the normalized % direction vectors to the first and to the i'thsatellites. % Last element is difference between 1st and i'th innerproducts % of direction vectors with satellite velocity vectors. for i =1:4, H(i,1:3) = DV(:,i+1)′ − DV(:,1)′; H(i,4) = DV(:,1)′ * V(:,1) −DV(:,i+1)′ * V(:,i+1); end H4(:,loop) = H(:,4); % step 6 % compute G via1 iteration % G = G*(2*eye(4)−H*G); G = inv(H); G4(loop,:) = G(4,:); %step 7 % compute dU as G * r_diff_error(:,loop) dU4 = G *r_diff_error(:,loop); dU = dU4(1:3); dt = dU4(4); % step 8 % updateposition and time estimate Ub = Ub + dU, time = time + dtans_vec(:,loop) = [Ub*,time]′; Ub = Ub + dU; time = time + dt; loop; end% while U_covar = G * r_diff_covar * G′; U_sigma(1) =sqrt(U_covar(1,1)); U_sigma(2) = sqrt(U_covar(2,2)); U_sigma(3) =sqrt(U_covar(3,3)); U_sigma(4) = sqrt(U_covar(4,4)); % Plot x,y,z or tconvergence vs loop iteration number. % 1 −> x, 2 −> y, 3 −> z, 4 −>time. % | % | % v plot(ans_vec(1,:)); end function [vec_satellite_ecf] =satellite_ephemeris(input_vector) % The following code for calculatingsatellite positions from ephemeris % data is taken from ICD-GPS-2000,table 20-IV. t = input_vector(1,1); m_0 = input_vector(2,1); delta_n =input_vector(3,1); c = input_vector(4,1); sqrt_a = input_vector(5,1);OMEGA_0 = input_vector(6,1); i_0 = input_vector(7,1); omega =input_vector(8,1); OMEGA_dot = input_vector(9,1); i_dot =input_vector(10,1); c_uc = input_vector(11,1); c_us =input_vector(12,1); c_rc = input_vector(13,1); c_rs =input_vector(14,1); c_ic = input_vector(15,1); c_is =input_vector(16,1); t_oe = input_vector(17,1); iode =input_vector(18,1); % Earth mass times universal gravity constant (seeRef 1, pp 34). mu = 3.986005e14; % Earth rotation rate. OMEGA_e_dot =7.2921151467e−5; % Orbit ellipse semimajor axis length a =(sqrt_a){circumflex over ( )}2; % Mean angular motion (Average angularrate of orbit). n_0 = (mu/(a{circumflex over ( )}3)){circumflex over( )}0.5; % Time elapsed since ephemeris reference epoch. t_k = t − t_oe;if t_k > 302400 t_k = t_k−604800; % correct time for EOW crossoverselseif t_k <− 302400 t_k = t_k+604800; % correct time for EOW crossoversend % Corrected mean angular motion. n = n_0 + delta_n; % Mean anomaly(false sat orbit angle in orbit plane, using mean motion). m_k = m_0 +n*t_k; % Eccentric anomaly (true sat orbit angle in orbit plane, fromellipse center). e_k = kepler(e,m_k,0.01,10000); % Terms for findingtrue anomaly. cosf_k = (cos(e_k) − e)/(1 − e*cos(e_k)); sinf_k = ((1 −e{circumflex over ( )}2){circumflex over ( )}0.5)*sin(e_k)/(1 −e*cos(e_k)); % True anomaly (true sat orbit angle in orbit plane, fromearth center). f_k = atan2(sinf_k,cosf_k); % Argument of latitude (satorbit angle from ascending node point on equator) phi_k = f_k + omega;sin2phik = sin(2*phi_k); cos2phik = cos(2*phi_k); % Argument of latitudecorrection. delta_u_k = c_us*sin2phik + c_uc*cos2phik; % Radiuscorrection. delta_r_k = c_rc*cos2phik + c_rs*sin2phik; % Correction toinclination. delta_i_k = c_ic*cos2phik + c_is*sin2phik; % Correctedargument of latitude. u_k = phi_k + delta_u_k; % Corrected radius. r_k =a*(1 − e*cos(e_k)) + delta_r_k; % Corrected inclination. i_k = i_0 +delta_i_k + i_dot*t_k; % Final x-y position in orbital plane. x_k_prime= r_k*cos(u_k); y_k_prime = r_k*sin(u_k); % Longitude of ascending node.omega_k = OMEGA_0 + (OMEGA_dot−OMEGA_e_dot)*t_k − OMEGA_e_dot*t_oe; %Final sat position in earth-fixed coordinates. x_k =x_k_prime*cos(omega_k) − y_k_prime*cos(i_k)*sin(omega_k); y_k =x_k_prime*sin(omega_k) + y_k_prime*cos(i_k)*cos(omega_k); z_k =y_k_prime*sin(i_k); vec_satellite_ecf = [x_k;y_k;z_k]; end function x =kepler(a,b,epsilon,N) % KEPLER KEPLER(a,b,epsilon,N) returns theiterated solution to % x=a*sin(x)+b. Iterates until the new value iswithin epsilon % of the old value or until the number of iterations %reaches N. old_x = 0; new_x = b; iterations=0; while((abs(new_x−old_x)>epsilon) & (iterations<N)) old_x = new_x;%disp(new_x); new_x = a*sin(old_x)+b; iterations = iterations+1; end ifiterations<N x=new_x; %disp(‘Iterations:’); %disp(iterations); %disp(a);%disp(b); %disp(x); elseif iterations= =N disp(‘Fixed Point IterationFailed’); end function [param_vec] = sv7_eph( ) EphemerisSV07length =167;%bytes SV_PRN = 7; t = 483634.754; t_ephem = 483042; weeknum = 763;codeL2 = 1; L2Pdata = 0; SVacc_raw = 7; SV_health = 0; IODC = 120; T_GD= 1.39698e−09; t_oc = 485504; a_f2 = 0; a_f1 = 1.13687e−13; a_f0 =0.000697501; SVacc = 32; IODE = 120; fit_intv1 = 0; C_rs = −17.2812;delta_n = 4.39733e−09; M_0 = −1.6207; C_uc = −9.27597e−07; ecc =0.00659284; C_us = 9.46783e−06; sqrt_A = 5153.78; t_oe = 485504; C_ic =1.47149e−07; OMEGA_0 = 0.793991; C_is = −9.31323e−09; i_0 = 0.962958;C_rc = 197.938; omega = −2.69189; OMEGADOT = −7.93212e−09; IDOT =7.82175e−11; Axis = 2.65614e+07; n = 0.00014585; r1me2 = 0.999978;OMEGA_n = −34.6095; ODOT_n = −7.29291e−05; param_vec = [t; M_0; delta_n;ecc; sqrt_A; OMEGA_0; i_0; omega; OMEGADOT; IDOT; C_uc; C_us; C_rc;C_rs; C_ic; C_is; t_oe; IODE]; function [param_vec] = sv2_eph( )EphemerisSV02length = 167;%bytes SV_PRN = 2; t = 483634.754; t_ephem =483042; weeknum = 763; codeL2 = 1; L2Pdata = 0; SVacc_raw = 7; SV_health= 0; IODC = 165; T_GD = 9.31323e−10; t_oc = 489600; a_f2 = 0; a_f1 =−2.16005e−12; a_f0 = −0.000109646; SVacc = 32; IODE = 165; fit_intv1 =0; C_rs = −30.6562; delta_n = 4.92235e−09; M_0 = 0.0785892; C_uc =−1.52551e−06; ecc = 0.0134761; C_us = 4.63426e−06; sqrt_A = 5153.66;t_oe = 489600; C_ic = 1.09896e−07; OMEGA_0 = −0.26998; C_is =1.63913e−07; i_0 = 0.953377; C_rc = 285.688; omega = −2.63668; OMEGADOT= −8.54964e−09; IDOT = −2.16795e−10; Axis = 2.65602e+07; n = 0.00014586;r1me2 = 0.999909; OMEGA_n = −35.9722; ODOT_n = −7.29297e−05; param_vec= [t, M_0; delta_n; ecc; sqrt_A; OMEGA_0; i_0; omega; OMEGADOT; IDOT;C_uc; C_us; C_rc; C_rs; C_ic; C_is; t_oe; IODE]; function [param_vec] =sv9_eph( ) EphemerisSV09length = 167;%bytes SV_PRN = 9; t = 483634.754;t_ephem = 483042; weeknum = 763; codeL2 = 1; L2Pdata = 0; SVacc_raw = 7;SV_health = 0; IODC = 592; T_GD = 1.39698e−09; t_oc = 489600; a_f2 = 0;a_f1 = −1.13687e−12; a_f0 = −5.16325e−06; SVacc = 32; IODE = 80;fit_intv1 = 0; C_rs = 80.1875; delta_n = 4.85127e−09; M_0 = −2.65977;C_uc = 4.25801e−06; ecc = 0.00292443; C_us = 5.34952e−06; sqrt_A =5153.73; t_oe = 489600; C_ic = 5.58794e−09; OMEGA_0 = −1.28659; C_is =2.23517e−08; i_0 = 0.950485; C_rc = 270.656; omega = −0.581175; OMEGADOT= −8.33535e−09; IDOT = 2.4501e−10; Axis = 2.6561e+07; n = 0.000145854;r1me2 = 0.999996; OMEGA_n = −36.9888; ODOT_n = −7.29295e−05; param_vec= [t; M_0; delta_n; ecc; sqrt_A; OMEGA_0; i_0; omega; OMEGADOT; IDOT;C_uc; C_us; C_rc; C_rs; C_ic; C_is; t_oe; IODE]; function [param_vec] =sv4_eph( ) EphemerisSV04length = 167;%bytes SV_PRN = 4; t = 483634.754;t_ephem = 483048; weeknum = 763; codeL2 = 1; L2Pdata = 0; SVacc_raw = 7;SV_health = 0; IODC = 711; T_GD = 1.39698e−09; t_oc = 489600; a_f2 = 0;a_f1 = 1.59162e−12; a_f0 = 4.03658e−05; SVacc = 32; IODE = 199;fit_intv1 = 0; C_rs = 15.4375; delta_n = 4.67448e−09; M_0 = 2.57307;C_uc = 9.76026e−07; ecc = 0.00303051; C_us = 6.10203e−06; sqrt_A =5153.62; t_oe = 489600; C_ic = 4.09782e−08; OMEGA_0 = 1.86007; C_is =1.11759e−08; i_0 = 0.963886; C_rc = 261.562; omega = −1.24408; OMEGADOT= −8.10534e−09; IDOT = 2.62511e−10; Axis = 2.65598e+07; n = 0.000145863;r1me2 = 0.999995; OMEGA_n = −33.8421; ODOT_n = −7.29293e−05; param_vec= [t; M_0; delta_n; ecc; sqrt_A; OMEGA_0; i_0; omega; OMEGADOT; IDOT;C_uc; C_us; C_rc; C_rs; C_ic; C_is; t_oe; IODE]; function [param_vec] =sv12_eph( ) % This ephemeris data was extracted from test0826.asc % byGlen Brooksby. EphemerisSV12length = 167;%bytes SV_PRN = 12; t =483634.754; t_ephem = 482748; weeknum = 763; codeL2 = 1; L2Pdata = 0;SVacc_raw = 2; SV_health = 0; IODC = 323; T_GD = 2.79397e−09; t_oc =489600; a_f2 = 0; a_f1 = −1.21645e−11; a_f0 = −2.44565e−05; SVacc = 4;IODE = 67; fit_intv1 = 0; C_rs = 42.75; delta_n = 1.96937e−09; m_0 =3.01393; C_uc = 2.28174e−06; ecc = 0.0147704; C_us = 5.3402e−06; sqrt_a= 5153.5; t_oe = 489600; C_ic = −1.65775e−07; omega_0 = −0.925903; C_is= 3.91155e−08; i_0 = 1.08737; C_rc = 338.719; omega = −0.167314;OMEGADOT = −6.63885e−09; IDOT = 1.2572e−10; Axis = 2.65585e+07; n =0.000145871; r1me2 = 0.999891; OMEGA_n = −36.6281; ODOT_n =−7.29278e−05; param_vec = [t; M_0; delta_n; ecc; sqrt_a; OMEGA_0; i_0;omega; OMEGADOT; IDOT; C_uc; C_us; C_rc; C_rs; C_ic; C_is; t_oe; IODE];

1. A method for identifying location of an object to be tracked,comprising: measuring data related to propagation time differencesbetween signals transmitted from a plurality of GPS satellites andreceived at said object to be tracked, said data comprising code wordphase measurements μ_(i) for a satellite at a time t_(R), whereμ_(i)=γ_(i)/T_(i) ^(C), and defined as time elapsed to time t_(R) fromthe beginning of a code word in the signal from satellite i in whicht_(R) falls, y_(i) being defined as the receiver code-time offset forsatellite i and T_(i) ^(C) being defined as the code period forsatellite i at time t_(R) in the signal received from satellite i, saidcode word phase measurements being simultaneously derived from thesignals transmitted from said plurality of satellites and received atthe object to be tracked; transmitting said data to a central station;and calculating at said central station the location of said object tobe tracked based upon the transmitted data and data derived from atleast one receiver apart from said object to be tracked receiving saidsignals from said plurality of satellites.
 2. The method of claim 1wherein the data transmitted to the central station includes satelliteidentification data so that the step of calculating the location of saidobject to be tracked is thereupon based further upon the satelliteidentification data.
 3. The method of claim 1 wherein the step ofcalculating the location of said object comprises calculating a point ofintersection of curves defined by said propagation time differences. 4.The method of claim 1 including the step of transmitting time signals tosaid object to be tracked over a separate channel so as to maintainclock accuracy at said object to be tracked.
 5. A method for identifyinglocation of an object to be tracked comprising: measuring data relatedto propagation time differences between signals transmitted from aplurality of GPS satellites and received at said object to be tracked,said data comprising bit phase measurements μ_(i) for a satellite i at atime t_(R), where μ_(i)=β_(i)/T_(i) ^(B), β_(i) being the receiverbit-time offset for satellite i and defined as time elapsed to timet_(R) from the beginning of a code word in the signal from satellite iin which t_(R) falls, T_(i) ^(B) being defined as the bit period forsatellite i at time t_(R) in the signal received from satellite i, saidbit phase measurements being simultaneously derived from the signalstransmitted from said plurality of satellites and received at the objectto be tracked; transmitting said data to a central station; andcalculating at said central station the location of said object to betracked based upon the transmitted data and data derived from at leastone receiver apart from said object to be tracked receiving said signalsfrom said plurality of satellites.
 6. The method of claim 5 and furtherincluding: recording, at said object, the time at which the data aresimultaneously derived; and transmitting the recorded time to saidcentral station.
 7. The method of claim 5 and further including:measuring, at said object to be tracked, delay between the time at whichthe data are recorded and the time when the data are transmitted to thecentral station; and transmitting the measured delay to said centralstation.
 8. The method of claim 5 wherein the data transmitted to thecentral station includes satellite identification data and the step ofcalculating the location of said object to be tracked is additionallybased upon the satellite identification data.
 9. The method of claim 5wherein the step of calculating the location of said object comprisescalculating a point of intersection of curves defined by saidpropagation time differences.
 10. The method of claim 5 wherein thesignals from said GPS satellites are received at said object to betracked, and including the step of transmitting time signals to saidobject to be tracked over a separate channel so as to maintain clockaccuracy at said object to be tracked.
 11. A method for identifyinglocation of an object to be tracked comprising: measuring data relatedto propagation time differences between signals transmitted from aplurality of satellites and received at said object to be tracked, saiddata comprising receiver code-time offsets for a satellite i and definedas time elapsed to a time t_(R) from the beginning of a code word in thesignal from satellite i in which t_(R) falls, and code periods in thesignal received from satellite i in which time t_(R) falls, saidplurality of satellites comprising GPS satellites, and including theadditional step of simultaneously deriving said receiver code-timeoffsets and code periods from signals received from the plurality ofsatellites at said object to be tracked; transmitting said data to acentral station; and calculating at said central station the location ofsaid object to be tracked based upon the transmitted data and dataderived from at least one receiver apart from said object to be trackedreceiving said signals from said plurality of satellites.
 12. A methodfor identifying location of an object to be tracked, comprising:measuring data related to propagation time differences between signalstransmitted from at least four GPS satellites and received at saidobject to be tracked, said data related to propagation time differencescomprising bit phase measurements simultaneously derived from saidsignals; transmitting said data, including satellite identificationdata, to a central station; measuring, at said object to be tracked,delay between the time at which the data are recorded and the time whenthe data are transmitted to the central station; transmitting themeasured delay to said central station; and calculating at said centralstation the location of said object to be tracked based upon thetransmitted data, the satellite identification data, and data derivedfrom at least one receiver apart from said object to be trackedreceiving said signals from said plurality of satellites, thecalculating step comprising: assuming a feasible value for acommunication time delay required for a signal transmitted from saidobject to be tracked to reach the central station; calculating thelocation of said object to be tracked based upon the satelliteidentification data and the assumed value of said communication timedelay; calculating a new value for said communication time delay basedupon the calculated location of said object to be tracked; andcalculating a corrected location of said object to be tracked based uponthe calculated new value for said communication time delay.
 13. Themethod of claim 12 including the additional step of: iterativelyrepeating the calculating steps until little change in location of saidobject to be tracked is observable.
 14. A system for identifyinglocation of an object to be tracked, comprising: means for measuringdata related to propagation time differences between signals transmittedfrom a plurality of GPS satellites and received at said object to betracked, each of said signals identifying an associated satellite, saidobject to be tracked including: receiver means for receiving signalsfrom at least four GPS satellites; and first processor means forprocessing data from the receiver means at predetermined time intervalsin synchronism with received signal events, said data being related topropagation time differences for said signals; receiver means apart fromsaid object for receiving said signals transmitted from said pluralityof satellites; a central station for calculating the location of saidobject based upon the measured data, data derived from said receivermeans apart from said object, and the satellite identification data; andtransmission means for transmitting the processed data to said centralstation; said system further including: second processor means at saidcentral station for determining the location of said object based on thedata received from said transmission means and data derived from saidreceiver means apart from said object.
 15. The system of claim 14wherein said signal events comprise a telemetry-word preamble signalevent in a GPS signal.
 16. The system of claim 15 wherein said firstprocessor means further comprises means for decoding a satellite timestamp from a predetermined one of the received GPS signals, based uponthe telemetry-word preamble signal event.
 17. A method system foridentifying location of an object to be tracked comprising: means formeasuring data related to propagation time differences between signalstransmitted from a plurality of GPS satellites and received at saidobject to be tracked, each of said signals identifying an associatedsatellite, said object to be tracked including: receiver means forreceiving signals from at least four GPS satellites; and first processormeans for calculating a receiver bit phase for each of said satellites,said bit phase for any satellite i at a time t_(R) being defined asβ_(i)/T_(i) ^(B), β_(i) being the receiver bit-time offset for satellitei and defined as time elapsed to time t_(R) from the beginning of a codeword in the signal from satellite i in which t_(R) falls, and T_(i) ^(B)being defined as the bit period for satellite i at time T₈ in the signalreceived from satellite i; receiver means apart from said object forreceiving said signals from said plurality of satellites; a centralstation; and transmission means for transmitting the calculated bitphases to said central station; said system further including: secondprocessor means at said central station for determining signalpropagation times between said plurality of satellites and said objectand for determining location of said object based upon the bit phasestransmitted by said transmission means and data derived from saidreceiver means apart from said object.
 18. A system for identifyinglocation of an object to be tracked, comprising: means for measuringdata related to propagation time differences between signals transmittedfrom a plurality of GPS satellites and received at said object to betracked, each of said signals identifying an associated satellite, saidobject to be tracked including: receiver means for receiving signalsfrom at least four GPS satellites, and first processor means forcalculating a bit-time offset for each of said satellites and fordetermining a bit period for each signal received from said satellites,said bit-time offset for a satellite i being defined as time elapsed toa time t_(R) from the beginning of a code word in the signal fromsatellite i in which t_(R) falls, said bit period for satellite i beingdetermined at time t_(R) in the signal from satellite i; receiver meansapart from said object for receiving said signals from said plurality ofsatellites; a central station; and transmission means for transmittingtime stamps, the calculated bit-time offsets and bit periods, andsatellite identification data, to said central station; said systemfurther including: second processor means at said central station fordetermining signal propagation times between said plurality ofsatellites and said object and for determining location of said objectbased upon the bit-time offsets and periods, time stamps, satelliteidentification data transmitted by said transmission means, and dataderived from said receiver means apart from said object.
 19. A systemfor identifying location of an object to be tracked comprising: areceiver for measuring data related to propagation time differencesbetween signal transmitted from a plurality of satellites and receivedat said object to be tracked, said data comprising receiver code-timeoffsets for a satellite i and defined as time elapsed to a time t _(R)from the beginning of a code word in the signal from satellite i inwhich t _(R) falls, and code periods in the signal received fromsatellite i in which time t _(R) falls, said plurality of satellitescomprising GPS satellites, said receiver further being adapted tosimultaneously derive said receiver code-time offsets and code periodsfrom signals received from the plurality of satellites at said object tobe tracked; and processor apparatus for storing the measured data for apredetermined length of time and for calculating, from the measureddata, the location of said object to be tracked; said receiver beingadapted to be deactivated while the measured data undergo processing.20. The system of claim 19 wherein the processor apparatus comprises afirst processor for storing the measured data for said predeterminedlength of time, and a second processor for calculating, from the storedmeasured data, the location of said object to be tracked; wherein thefirst processor has a means for processing data from the receiver meansat predetermined time intervals in synchronism with received signalevents, the data being related to propagation time differences for saidsignal.